{ Image Macros for intensity calculations on Photometrics CCD images }

{Globals}
var
    {16 to 8 conversion info}
    ymin, ymax,
    
    {pid numbers for images}
    customLUTpid, smoothPid,
    raw16Pid, rawWidth, rawHeight,
    proc16Pid, proc8Pid, procXmin, procXmax,
    temp16Pid, flagPid, 
    mask1Pid, maskWidth, maskHeight,
    seg8aPid, seg8bPid, bkgSegPid,
    dark16Pid, unif16Pid,
    stackPid,

    {segmentation}
    segN, molecN,
    
    {YES/NO parameters}
    skipFlat, {nonzero = skip flat field correction}
    skipMedian, {nonzero = skip median filter step}
    skipSmooth, {nonzero = skip smooth step}

    {smooth kernel parameters}
    {flat field ratio constant}
    
: integer;
{initialize/restore globals}
begin
    requiresUser('Pixel16u',2);
    requiresUser('GetPutPixel',1);
    requiresUser('timer',1);
    requiresUser('utilities',1);
    requiresUser('markup',1);
    
    ymin := getMemo('ymin');
    ymax := getMemo('ymax');

    {pid numbers}

    customLUTpid := getMemo('customLUTpid');
    smoothPid := getMemo('smoothPid');
    raw16Pid := getMemo('raw16Pid');
    rawWidth := getMemo('rawWidth');
    rawHeight := getMemo('rawHeight');
    proc16Pid := getMemo('proc16Pid');
    temp16Pid := getMemo('temp16Pid');
    flagPid := getMemo('flagPid');
    mask1Pid := getMemo('mask1Pid');
    maskWidth := getMemo('maskWidth');
    maskHeight := getMemo('maskHeight');
    proc8Pid := getMemo('proc8Pid');
    procXmin := getMemo('procXmin');
    procXmax := getMemo('procXmax');
    seg8aPid := getMemo('seg8aPid');
    seg8bPid := getMemo('seg8bPid');
    bkgSegPid := getMemo('bkgSegPid');
    dark16Pid := getMemo('dark16Pid');
    unif16Pid := getMemo('unif16Pid');
    stackPid := getMemo('stackPid');
    
    segN := getMemo('segN');
    if segN < 1 then segN := 1;
    if segN > 250 then segN := 250;
    setMemo('segN',segN);

    molecN := getMemo('molecN');
    if molecN < 1 then molecN := 1;
    if molecN > 250 then molecN := 250;
    setMemo('molecN',molecN);

    SetBackgroundColor(0);
    SetForeGroundColor(255);
end;

procedure disposePic(pid: integer);
begin
    selectPic(pid);
    dispose;
end;

procedure autoDispose(p);
begin
    if pidExists(p) then begin
        disposePic(p);
    end;
end;
procedure createSmoothKernel;
var
    x, y: integer;
    sum: real;
begin
    RequiresUser('getputpixel', 1);
    AutoDispose(smoothPid);
    SaveState;
    SetNewSize(kw * 4, kh + 1);
    MakeNewWindow('kernel ', kx : 0, 'x ', ky : 0, 'y ', kw : 0, 'w ', kh : 0, 'h');
    smoothPid := pidNumber;
    SetMemo('smoothPid', smoothPid);
    RestoreState;
    putPixel(0, 0, kx);
    putPixel(1, 0, ky);
    putPixel(2, 0, kw);
    putPixel(3, 0, kh);
    sum := 0.0;
    for x := -kx to kw - kx - 1 do begin
        for y := -ky to kh - ky - 1 do begin
            sum := sum + exp(-sqrt(sqr(x) + sqr(y)));
        end;
    end;
    sum := 32000.0 / sum; {nearly maximum before overflow on 65535 pixel}
    for x := -kx to kw - kx - 1 do begin
        for y := -ky to kh - ky - 1 do begin
            putPixel32s(smoothPid, x + kx, y + ky + 1, sum * exp(-sqrt(sqr(x) + sqr(y))));
        end;
    end;
    sum := 0.0;
    for x := -kx to kw - kx - 1 do begin
        for y := -ky to kh - ky - 1 do begin
            sum := sum + getPixel32s(smoothPid, x + kx, y + ky + 1);
        end;
    end;
    showmessage('kernel sum = ', sum);
    SelectPic(smoothPid);
    MakeRoi(0, 1, kw * 4, kh);
end;

procedure stdSmooth;
var
    kx, ky, kw, kh: integer;
begin
    kx := 1;
    ky := 1;
    kw := 3;
    kh := 3;
    createSmoothKernel;
end;

macro 'Create standard 3x3 smoothing kernel';
begin
    stdSmooth;
end;

macro 'Create arbitrary smoothing kernel';
var
    kx, ky, kw, kh: integer;
begin
    kw := GetNumber('kernel width',5);
    kh := GetNumber('kernel height',5);
    kx := GetNumber('kernel x center',kw div 2);
    ky := GetNumber('kernel y center',kh div 2);
    createSmoothKernel;
end;

procedure setMaskSize;
var
    kx, ky, kw, kh: integer;
    front: integer;
begin
    front := pidNumber;
    {find size of smoothing kernel}
    ChoosePic(smoothPid);
    kx := getPixel(0, 0);
    ky := getPixel(1, 0);
    kw := getPixel(2, 0);
    kh := getPixel(3, 0);

    MakeRoi(0, 1, kw * 4, kh);
    maskWidth := rawWidth + kw - 1;
    maskHeight := rawHeight + kh - 1;
    SetMemo('maskWidth',maskWidth);
    SetMemo('maskHeight',maskHeight);
    ChoosePic(front);
end;

procedure forceROIWithin;
var
    left, top, rwidth, rheight, iwidth, iheight: integer;
begin
    GetPicSize(iwidth,iheight);
    GetRoi(left,top,rwidth,rheight);
    if rwidth = 0 then selectAll; {this fixes most cases}
    GetRoi(left,top,rwidth,rheight);
    if (left < 0) 
    or (top < 0) 
    or (left + rwidth > iwidth) 
    or (top + rheight > iheight) then begin
        putmessage('ROI must not extend outside image');
        exit;
        {make ROI doesn't hack it if ROI wasn't rectangular...}
        if left < 0 then begin
            rwidth := rwidth + left;
            left := 0;
        end;
        if top < 0 then begin
            rheight := rheight + top;
            top := 0;
        end;
        if left + rwidth > iwidth then begin
            rwidth := iwidth - left;
        end;
        if top + rheight > iheight then begin
            rheight := iheight - top;
        end;
        makeroi(left,top,rwidth,rheight);
    end;
end;

procedure forceUncalib;
begin
    SelectPic(proc8Pid);
    if Calibrated then begin
        selectAll;
        copy;
        disposePic(proc8Pid);
        SaveState;
        setNewSize(rawWidth, rawHeight);
        makeNewWindow('Processed 8 bit image');
        proc8Pid := pidNumber;
        setMemo('proc8Pid',proc8Pid);
        RestoreState;
        Paste;
        KillRoi;
    end;
    SelectPic(proc8Pid);
end;

{adjust xmin/xmax using mean ± stdev}
procedure enhanceStdev;
var
    mean, sigma, coef: real;
begin
    choosePic(proc16Pid);
    KillROI;
    choosePic(proc8Pid);
    forceROIWithin;
    forceUncalib;
    KillROI;
    coef := (procXmax - procXmin + 1) / (ymax - ymin + 1);
{might not work if coef < 0???}
    linLUT16uto8(customLUTPid, procXmin, procXmax, ymin, ymax);
    Cnvrt16uto8(proc16Pid, customLUTPid, proc8Pid);
    RestoreRoi; {take mean & stdev over ROI of 8 bit image}
    SetOptions('Area,Mean,Std. Dev.,User1,User2');
    SetCounter(513);
    Measure;
    mean := (rmean[rCount]-ymin) * coef + procXmin + coef / 2;
    sigma := rStdDev[rCount] * coef + coef / 2;
    ruser1[rCount] := coef;
    ruser2[rCount] := mean;
    {serious round off errors happen when sigma < coef }
    {so that the mean is not known well enough, }
    {image comes out white or black}
    if sigma < coef then sigma := coef;
    procXmin := mean - 2*sigma; {this needs to be an adjustable parameter}
    procXmax := mean + 4*sigma;
    SetMemo('procXmin',procXmin);
    SetMemo('procXmax',procXmax);
end;

{adjust xmin/xmax using mean and maximum}
procedure adjMeanMax;
var
    mean, sigma, coef: real;
begin
    choosePic(proc16Pid);
    KillROI;
    choosePic(proc8Pid);
    forceROIWithin;
    forceUncalib;
    KillROI;
    coef := (procXmax - procXmin + 1) / (ymax - ymin + 1);
{might not work if coef < 0???}
    linLUT16uto8(customLUTPid, procXmin, procXmax, ymin, ymax);
    Cnvrt16uto8(proc16Pid, customLUTPid, proc8Pid);
    RestoreRoi; {take mean & stdev over ROI of 8 bit image}
    SetCounter(513);
    SetOptions('Area,Mean,Std. Dev.,User1,User2');
    Measure;
    InsetRoi(-20);
    Measure;
    mean := (rmean[rCount]-ymin) * coef + procXmin + coef / 2;
{sigma is not standard deviation but rather adjusted to put mean at 1/4 and max at max}
    sigma := (rmax[rCount-1]-ymin) * coef + procXmin + coef / 2;
    sigma := (sigma - mean) / 4;
    ruser1[rCount] := coef;
    ruser2[rCount] := mean;
    if sigma < coef then sigma := coef;
    procXmin := mean - 0.75*sigma; {this needs to be an adjustable parameter}
    procXmax := mean + 4*sigma;
    SetMemo('procXmin',procXmin);
    SetMemo('procXmax',procXmax);
end;


{display 16 bit data into the 8 bit image using specified xmin/xmax}
procedure show16;
var
    lower, upper: integer;
begin
    choosePic(proc16Pid);
    KillROI;
    choosePic(proc8Pid);
    forceROIWithin;
    KillROI;
    linLUT16uto8(customLUTPid, procXmin, procXmax, ymin, ymax);
    Cnvrt16uto8(proc16Pid, customLUTPid, proc8Pid);
    RestoreRoi;
    SelectPic(proc8Pid);
    GetThresholds(lower, upper);
    ShowMessage(procXmin,' min\',procXmax,' max\',lower,' lower\',upper,' upper\');
end;

macro 'Set 8 bit display range';
begin
    ymin := getnumber('gray level for smallest pixel value',ymin);
    ymax := getnumber('gray level for largest pixel value',ymax);
    SetMemo('ymin',ymin);
    SetMemo('ymax',ymax);
end;

procedure swapTemp16;
var
    temp: integer;
begin
    temp := temp16Pid;
    temp16Pid := proc16Pid;
    proc16Pid := temp;
    SetMemo('temp16Pid',temp16Pid);
    SetMemo('proc16Pid',proc16Pid);
    temp := pidNumber;
    choosePic(proc16Pid);
    SetPicName('Processed 16 bit image');
    killROI;
    choosePic(temp16Pid);
    SetPicName('Temporary 16 bit image');
    killROI;
    choosePic(temp);
end;

procedure hide8image;
var
    width, height: integer;
begin
    SelectPic(proc8Pid);
    setDensitySlice(0,0);
    setforegroundcolor(255);
    setbackgroundcolor(0);
    selectAll;
    clear;
    getPicSize(width, height);
    moveto(width div 3, height div 3);
    writeln('Press 8 to display image');
end;

procedure checkSize(p,w,h: integer);
var
    width, height, front: integer;
begin
    if pidExists(p) then begin
        front := pidNumber;
        choosePic(p);
        getPicSize(width, height);
        choosePic(front);
        if (width <> w) or (height <> h) then disposePic(p);
    end;
end;

{if the scratch windows are wrong size or missing, create them}
procedure makeScratchIfNeed;
var
    width, height: integer;
begin
    SelectPic(pidNumber);
    saveState;
    if (ymin = 0) and (ymax = 0) then begin
        ymin := 1;
        ymax := 254;
    end;
    if (ymin < 0) or (ymin > 255) then ymin := 1;
    if (ymax < 0) or (ymax > 255) then ymax := 254;
    if ymin > ymax then begin
        ymin := 1;
        ymax := 254;
    end;
    SetMemo('ymin',ymin);
    SetMemo('ymax',ymax);
    if not pidExists(customlutPid) then begin
        setNewSize(256,256);
        makeNewWindow('custom LUT');
        SelectAll;
        KillRoi;
        customLUTpid := pidNumber;
        SetMemo('customLUTpid',customLUTpid);
    end;
    linLUT16uto8(customLUTpid, 0, 65535, ymin, ymax);

    if not pidExists(smoothPid) then begin
        stdSmooth;
    end;

    if not pidExists(raw16Pid) then begin
        putMessage('makeScratch no raw16');
        exit;
    end;
    choosePic(raw16Pid);
    getPicSize(width, height);
    rawWidth := (width div 4) * 2;
    rawHeight := height;
    setMemo('rawWidth',rawWidth);
    setMemo('rawHeight',rawHeight);
    if rawWidth * 2 <> width then begin
        putMessage('makeScratch raw width not multiple of 4');
        exit;
    end;

    checkSize(proc16Pid,rawWidth * 2,rawHeight);
    if not pidExists(proc16Pid) then begin
        setNewSize(rawWidth * 2, rawHeight);
        makeNewWindow('Processed 16 bit image');
        SelectAll;
        KillRoi;
        proc16Pid := pidNumber;
        SetMemo('proc16Pid',proc16Pid);
    end;

    checkSize(temp16Pid,rawWidth * 2,rawHeight);
    if not pidExists(temp16Pid) then begin
        setNewSize(rawWidth * 2, rawHeight);
        makeNewWindow('Temporary 16 bit image');
        SelectAll;
        KillRoi;
        temp16Pid := pidNumber;
        SetMemo('temp16Pid',temp16Pid);
    end;
    
    checkSize(dark16Pid,rawWidth * 2,rawHeight);
    if not pidExists(dark16Pid) then begin
        setNewSize(rawWidth * 2, rawHeight);
        makeNewWindow('Dark 16 bit image');
        SelectAll;
        KillRoi;
        dark16Pid := pidNumber;
        SetMemo('dark16Pid',dark16Pid);
    end;
    
    checkSize(unif16Pid,rawWidth * 2,rawHeight);
    if not pidExists(unif16Pid) then begin
        setNewSize(rawWidth * 2, rawHeight);
        makeNewWindow('Uniform 16 bit image');
        SelectAll;
        KillRoi;
        unif16Pid := pidNumber;
        SetMemo('unif16Pid',unif16Pid);
    end;
    
    checkSize(flagPid,rawWidth,rawHeight);
    if not pidExists(flagPid) then begin
        setNewSize(rawWidth, rawHeight);
        makeNewWindow('smoothing flag image');
        SelectAll;
        KillRoi;
        flagPid := pidNumber;
        SetMemo('flagPid',flagPid);
    end;
    
    setMaskSize;
    checkSize(mask1Pid,maskWidth,maskHeight);
    if not pidExists(mask1Pid) then begin
        setNewSize(maskWidth, maskHeight);
        makeNewWindow('smoothing mask image 1');
        SelectAll;
        KillRoi;
        mask1Pid := pidNumber;
        SetMemo('mask1Pid',mask1Pid);
    end;

    checkSize(proc8Pid,rawWidth,rawHeight);
    if not pidExists(proc8Pid) then begin
        setNewSize(rawWidth, rawHeight);
        makeNewWindow('Processed 8 bit image');
        SelectAll;
        KillRoi;
        proc8Pid := pidNumber;
        setMemo('proc8Pid',proc8Pid);
    end;

    checkSize(seg8aPid,rawWidth,rawHeight);
    if not pidExists(seg8aPid) then begin
        setNewSize(rawWidth, rawHeight);
        makeNewWindow('Segments A');
        SelectAll;
        KillRoi;
        seg8aPid := pidNumber;
        setMemo('seg8aPid',seg8aPid);
    end;

    checkSize(seg8bPid,rawWidth,rawHeight);
    if not pidExists(seg8bPid) then begin
        setNewSize(rawWidth, rawHeight);
        makeNewWindow('Segments B');
        SelectAll;
        KillRoi;
        seg8bPid := pidNumber;
        setMemo('seg8bPid',seg8bPid);
    end;

    checkSize(bkgSegPid,rawWidth,rawHeight);
    if not pidExists(bkgSegPid) then begin
        setNewSize(rawWidth, rawHeight);
        makeNewWindow('Background Segments');
        SelectAll;
        KillRoi;
        bkgSegPid := pidNumber;
        setMemo('bkgSegPid',bkgSegPid);
    end;

    restoreState;
end;

procedure copyRawToProc;
begin
    choosePic(raw16Pid);
    selectAll;
    copy;
    killRoi;
    choosePic(proc16Pid);
    selectAll;
    paste;
    killRoi;
end;

{
import arbitrary IPLab image:
if (getpixel(0,0) <> ord('2')) 
or (getpixel(1,0) <> ord('.'))
or (getpixel(2,0) <> ord('3'))
or (getpixel(3,0) <> ord('a'))
or (getpixel(4,0) <> 0)
or (getpixel(5,0) <> 1) {short int}
then
else
 width := ((getpixel(6,0) * 256 + getpixel(7,0)) * 256 + getpixel(8,0)) * 256 + getpixel(9,0);
 height := ((getpixel(10,0) * 256 + getpixel(11,0)) * 256 + getpixel(12,0)) * 256 + getpixel(13,0);
 offset := 2120;
}

{if there is no raw data image, import one}
procedure importIfNeed;
var 
   origPid: integer;
begin
    origPid := 0;
    if not pidExists(raw16Pid) then begin
        SaveState;
        SetImport('8-bits,Custom');
        SetCustom(2634,1034,2124);
        Import('');
        origPid := pidNumber;
        {MakeNewWindow will not make odd width windows.}
        {Therefore, 16 bit images must be even # pixels wide}
        {or width multiple of 4}
        SetNewSize(2632,1032);
        MakeRoi(0, 2, 2632,1032);
        Copy;
        MakeNewWindow(GetPicName);
        raw16Pid := pidNumber;
        SetMemo('raw16Pid',raw16Pid);
        Paste;
        KillROI;
        disposePic(origPid);
        RestoreState;
    end;
    makeScratchIfNeed;
    SelectPic(proc8Pid);
    if origPid <> 0 then begin
        CopyRawToProc;
    end;
end;

macro '[1] copy proc to dark image';
begin
    SelectPic(proc8Pid);
    choosePic(proc16Pid);
    selectAll;
    copy;
    killRoi;
    choosePic(dark16Pid);
    selectAll;
    paste;
    killRoi;
    SelectPic(proc8Pid);
end;
macro '[2] copy proc to uniform image';
begin
    SelectPic(proc8Pid);
    choosePic(proc16Pid);
    selectAll;
    copy;
    killRoi;
    choosePic(unif16Pid);
    selectAll;
    paste;
    killRoi;
    SelectPic(proc8Pid);
end;

macro '[d] subtract dark image';
begin
    SelectPic(proc8Pid);
    choosePic(proc16Pid);
    killROI;
    choosePic(temp16Pid);
    killROI;
    choosePic(dark16Pid);
    killROI;
    sub16u(proc16Pid,dark16Pid,temp16Pid);
    swapTemp16;
    hide8Image;
    writeln('subtract dark');
end;

macro '[f] flat field -- divide by uniform image';
begin
    SelectPic(proc8Pid);
    choosePic(proc16Pid);
    killROI;
    choosePic(temp16Pid);
    killROI;
    choosePic(unif16Pid);
    killROI;
    ratio16u(proc16Pid,unif16Pid,temp16Pid,16384);
    swapTemp16;
    hide8Image;
    writeln('divide by uniform');
end;

{macro 'include median filter step';
macro 'skip median filter step';
macro 'include smoothing step';
macro 'skip smoothing step';
macro 'include flat field step';
macro 'skip flat field step';}

macro '[a] start over from raw image';
begin
    SelectPic(proc8Pid);
    importIfNeed;
    CopyRawToProc;
    hide8Image;
    writeln('raw data');
end;

macro '[z] undo last 16 bit transform';
begin
    SelectPic(proc8Pid);
    swapTemp16;
    hide8Image;
    writeln('undo');
end;

procedure doReduceNoise;
begin
    SelectPic(proc8Pid);
    {actually only need to copy the border}
    choosePic(proc16Pid);
    selectAll;
    copy;
    killRoi;
    choosePic(temp16Pid);
    selectAll;
    paste;
    {end copy}
    choosePic(proc16Pid);
    makeRoi(2,1,(rawWidth-2)*2,rawHeight-2);
    choosePic(temp16Pid);
    makeRoi(2,1,(rawWidth-2)*2,rawHeight-2);
    median16u(proc16Pid,temp16Pid);
    choosePic(proc16Pid);
    killROI;
    choosePic(temp16Pid);
    killROI;
    swapTemp16;
    hide8Image;
    writeln('reduce noise');
end;

macro '[r]reduce noise';
begin
    doReduceNoise;
end;

procedure doRadMed(radius: real);
var
    r: integer;
begin
    SelectPic(proc8Pid);
    r := round(radius + 0.5);
    {actually only need to copy the border}
    choosePic(proc16Pid);
    selectAll;
    copy;
    killRoi;
    choosePic(temp16Pid);
    selectAll;
    paste;
    killRoi;
    {end copy}
    choosePic(proc16Pid);
    makeRoi(2*r,r,(rawWidth-2*r)*2,rawHeight-2*r);
    choosePic(temp16Pid);
    makeRoi(2*r,r,(rawWidth-2*r)*2,rawHeight-2*r);
    radMedian16u(proc16Pid,temp16Pid,radius);
    choosePic(proc16Pid);
    killROI;
    choosePic(temp16Pid);
    killROI;
    swapTemp16;
    hide8Image;
    writeln('radial median filter');
end;

macro 'radial median filter';
var
    radius: real;
begin
    radius := getNumber('radius',10);
    doRadMed(radius);
end;

macro '[4] Radial median filter 30 10';
begin
    doRadMed(30);
    doRadMed(10);
end;

macro '[m]min spatial filter';
begin
    SelectPic(proc8Pid);
    choosePic(proc16Pid);
    makeRoi(2,1,(rawWidth-2)*2,rawHeight-2);
    choosePic(temp16Pid);
    makeRoi(2,1,(rawWidth-2)*2,rawHeight-2);
    minspat16u(proc16Pid,temp16Pid);
    choosePic(proc16Pid);
    killROI;
    choosePic(temp16Pid);
    killROI;
    swapTemp16;
    hide8Image;
    writeln('min spatial');
end;

macro '[s] smooth';
var
    kx, ky, kw, kh: integer;
begin
    SelectPic(proc8Pid);
    SetBackgroundColor(0);
    ChoosePic(proc16Pid);
    killRoi;
    ChoosePic(temp16Pid);
    killRoi;
    ChoosePic(smoothPid);
    kx := getPixel(0, 0);
    ky := getPixel(1, 0);
    kw := getPixel(2, 0);
    kh := getPixel(3, 0);
    MakeRoi(0, 1, kw * 4, kh);
    ChoosePic(flagPid);
    SelectAll;
    Clear;
    KillRoi;
    ChoosePic(mask1Pid);
    MakeRoi(kx, ky, rawWidth, rawHeight);
    Clear;
    SetForegroundColor(255);
    MakeRoi(0, 0, kx, rawHeight + kh);
    Fill;
    MakeRoi(kx + rawWidth, 0, kw - kx - 1, rawHeight + kh);
    Fill;
    MakeRoi(kx, 0, rawWidth, ky);
    Fill;
    MakeRoi(kx, ky + rawHeight, rawWidth, kh - ky - 1);
    Fill;
    {Mask image must have an ROI same size as image and}
    {with borders matching kernel, thus:}
    MakeRoi(kx, ky, rawWidth, rawHeight);
    Convolve16u(flagPid, proc16Pid, smoothPid, kx, ky, mask1Pid, temp16Pid);
    swapTemp16;
    hide8Image;
    writeln('smooth');
end;

macro 'Copy 8 bit image to seg8b and make binary';
var
    lower, upper: integer;
begin
    SelectPic(proc8Pid);
    GetThresholds(lower, upper);
    SelectAll;
    Copy;
    ChoosePic(seg8bPid);
    SelectAll;
    Paste;
    if upper = 255 then begin
        SetThreshold(lower);
    end else begin
        SetDensitySlice(lower,upper);
    end;
    MakeBinary;
    SelectPic(seg8bPid);
end


macro '[3]Load a new image';
begin
    disposePic(raw16Pid);
    importIfNeed;
    minmax16u(proc16Pid, procXmin, procXmax);
    SetMemo('procXmin',procXmin);
    SetMemo('procXmax',procXmax);
    show16;
    enhanceStdev;
    show16;
end;

macro 'Front 16 bit image is raw data';
begin
    raw16Pid := pidNumber;
    SetMemo('raw16Pid',raw16Pid);
    SelectPic(raw16Pid);
    makeScratchIfNeed;
    CopyRawToProc;
    hide8Image;
    writeln('raw data');
end;

macro '[7]Convert to 8 bit with mean max scaling';
begin
    importIfNeed;
    minmax16u(proc16Pid, procXmin, procXmax);
    SetMemo('procXmin',procXmin);
    SetMemo('procXmax',procXmax);
    show16;
    adjMeanMax;
    show16;
end;


macro '[*]Convert to 8 bit with min max scaling';
begin
    importIfNeed;
    minmax16u(proc16Pid, procXmin, procXmax);
    SetMemo('procXmin',procXmin);
    SetMemo('procXmax',procXmax);
    show16;
end;

macro '[8]Convert to 8 bit with mean ± stdev scaling';
begin
    importIfNeed;
    minmax16u(proc16Pid, procXmin, procXmax);
    SetMemo('procXmin',procXmin);
    SetMemo('procXmax',procXmax);
    show16;
    enhanceStdev;
    show16;
end;

macro '[¥]Enhance ROI of 8 bit image';
begin
    importIfNeed;
    enhanceStdev;
    show16;
end;

macro '[9]reduce xmin';
begin
    importIfNeed;
    procXmin := round(procXmin - 0.1*(procXmax - procXmin) - 1);
    if procXmin > procXmax then procXmax := procXmin + 1;
    SetMemo('procXmin',procXmin);
    show16;
end;

macro '[»]increase xmin';
begin
    importIfNeed;
    procXmin := round(procXmin + 0.1*(procXmax - procXmin) + 1);
    if procXmin > procXmax then procXmax := procXmin + 1;
    SetMemo('procXmin',procXmin);
    show16;
end;

macro '[0]reduce xmax';
begin
    importIfNeed;
    procXmax := round(procXmax - 0.1*(procXmax - procXmin) - 1);
    if procXmax < procXmin then procXmin := procXmax - 1;
    SetMemo('procXmin',procXmin);
    show16;
end;
macro '[¼]increase xmax';
begin
    importIfNeed;
    procXmax := round(procXmax + 0.1*(procXmax - procXmin) + 1);
    if procXmin > procXmax then procXmin := procXmax - 1;
    SetMemo('procXmin',procXmin);
    show16;
end;

macro '[g] Show processed 8 bit image';
begin
    SelectPic(proc8Pid);
end;

procedure adjustsegN(offset: integer);
var
    wrap: integer;
begin
    if offset < 0 then
        wrap := 250
    else
        wrap := 1;
    segN := segN + offset;
    if segN > 250 then segN := wrap;
    if segN < 1 then segN := wrap;
    setMemo('segN',segN);
end;

macro '[j] Hilight next segment';
var
    oldsegN: integer;
begin
    adjustsegN(0);
    oldsegN := segN;
    SelectPic(seg8aPid);
    killRoi;
    SetCounter(513);
    measure;
    SetCounter(rCount - 1);
    repeat
        adjustsegN(1);
    until (histogram[segN] <> 0) or (segN = oldsegN);
    setDensitySlice(segN,segN);
    ShowMessage(segN,' Segment\');
end;

procedure appendROI;
var 
    fg, lower,upper: integer;
begin
    fg := pidNumber;
    GetThreshold(lower,upper);
    SetDensitySlice(0,0);
    KillRoi;
    RestoreRoi;
    Clear;
    ChoosePic(seg8aPid);
    SetBackgroundColor(0);
    SetForegroundColor(segN);
    RestoreRoi;
    fill;
    SetCounter(513);
    Measure;
    rX[segN+256] := rX[rCount];
    rY[segN+256] := rY[rCount];
    rMin[segN+256] := segN;
    rMax[segN+256] := molecN;
    SetForegroundColor(255);
    SelectPic(fg);
    ShowMessage(segN,' Segment\',lower,'lower\',upper,'upper\');
    if upper = 255 then
         SetThreshold(lower)
    else
         SetDensitySlice(lower,upper);
end;

macro '[n]ROI is next segment';
begin
    adjustsegN(1);
    appendROI;
end;

macro '[v]Next seg will be another molecule';
begin
    molecN := molecN + 1;
end;

procedure eraseSegment;
var 
    fg, lower,upper: integer;
begin
    fg := pidNumber;
    GetThreshold(lower,upper);
    SetDensitySlice(0,0);
    ChoosePic(seg8aPid);
    ChangeValues(segN,segN,0);
    rX[segN+256] :=0;
    rY[segN+256] :=0;
    rUser1[segN+256] := 0;
    rUser2[segN+256] := 0;
    rMin[segN+256] := 0;
    rMax[segN+256] := 0;
    SelectPic(fg);
    ShowMessage(segN,' Segment\',lower,'lower\',upper,'upper\');
    if upper = 255 then
         SetThreshold(lower)
    else
         SetDensitySlice(lower,upper);
end;

macro '[b] erase this segment';
begin
    eraseSegment;
end;

macro '[F1]Clear segments A and reset segment number';
var
    i: integer;
begin
    SelectPic(seg8aPid);
    SetBackgroundColor(0);
    SelectAll;
    Clear;
    SelectPic(proc8Pid);
    segN := 1;
    SetMemo('segN',segN);
    molecN := 1;
    SetMemo('molecN',molecN);
    for i := 256 to 511 do begin
        rX[i] := 0;
        rY[i] := 0;
        rUser1[i] := 0;
        rUser2[i] := 0;
        rMin[i] := 0;
        rMax[i] := 0;
        rAngle[i] := 0;
    end;
    SetCounter(513);
end;

macro '[F2]dilate segments A onto segments B';
var
    r: integer;
    width,height: integer;
begin
    SelectPic(seg8aPid);
    r := GetNumber('dilation radius',5);
    SelectAll;
    Copy;
    InsetRoi(r+1);
    choosePic(seg8bPid);
    SelectAll;
    Paste;
    InsetRoi(r+1);
    Dilate8Circular(seg8aPid, seg8bPid, r);
    choosePic(seg8aPid);
    KillRoi;
    SelectPic(seg8bPid);
    KillRoi;
end;

macro '[F3]copy segments B onto segments A';
begin
    SelectPic(seg8bPid);
    SelectAll;
    Copy;
    KillRoi;
    SelectPic(seg8aPid);
    SelectAll;
    Paste;
    KillRoi;
end;

macro '[F4]copy 8 bit image to Background Segments';
begin
    SelectPic(Proc8Pid);
    SelectAll;
    Copy;
    KillRoi;
    SelectPic(BkgSegPid);
    SelectAll;
    Paste;
    KillRoi;
    SetCounter(513);
    Measure;
    SetDensitySlice(rMean[rCount]-rStdDev[rCount]/2, rMean[rCount]+rStdDev[rCount]/2);
end;

macro '[F5]Tabulate intensities using segments A';
var
    i, j, k, r, x, area, sum, lastMol: integer;
    sumPid: integer;
    avgBkg, total: real;
begin
    showmessage('finding neighborhood of segments');
    SelectPic(seg8aPid);
    SelectAll;
    Copy;
    KillRoi;

    r := 5;
    ChoosePic(flagPid);
    SelectAll;
    Paste;
    ChangeValues(1,255,1);
    Copy;
    InsetRoi(r+1);
    ChoosePic(seg8bPid);
    SelectAll;
    Paste;
    InsetRoi(r+1);
    Dilate8Circular(flagPid, seg8bPid, r);
    ChoosePic(seg8bPid);
    SelectAll;
    Copy;
    r := 10;
    ChoosePic(flagPid);
    SelectAll;
    Paste;
    InsetRoi(r+1);
    ChoosePic(seg8bPid);
    SelectAll;
    Invert;
    Copy;
    Invert;
    InsetRoi(r+1);
    Dilate8Circular(seg8bPid, flagPid, r);
    ChoosePic(flagPid);
    SelectAll;
    Paste;
    DoAnd;
    ChangeValues(1,255,255);
    Showmessage('finding background pixels');
    ChoosePic(BkgSegPid);
    SelectAll;
    Copy;
    ChoosePic(flagPid);
    SelectAll;
    Paste;
    DoAnd;
    SelectAll;
    Copy;
    ChoosePic(Seg8bPid);
    Paste;
    ChoosePic(Seg8aPid);
    SelectAll;
    Copy;
    ChoosePic(Seg8bPid);
    Paste;
    DoOr;
    ShowMessage('calculating');

    SaveState;
    SetNewSize(32,32);
    MakeNewWindow('Sums');
    sumPid := pidNumber;
    
    Sum16uMark(proc16Pid,Seg8bPid,sumPid);
    {also need standard deviation, min, max}
    SetCounter(513);
    SetOptions('area,user1,user2,Min/Max,X-Y Center,Angle');
    ChoosePic(Seg8bPid);
    KillRoi;
    Measure; {get area from histogram array}
    SetCounter(513);
    i := 0;
    if histogram[255] = 0 then begin
        putMessage('no background');
        exit;
    end;
    avgBkg := GetPixVec32s(sumPid,255) / histogram[255];
    total := 0;
    lastMol := 0;
    rMax[511] := -1; {force normalize last molecule}
    for k := 1 to 255 do begin
        sum := GetPixVec32s(sumPid,k);
        area := histogram[k];
        if area <> 0 then begin
            if  rMax[k+256] <> lastMol then begin
                if total <> 0 then begin
                    for j := 1 to i do begin
                        if rMax[j] = lastMol then begin
                            rUser2[j] := rAngle[j] / total;
                        end;
                    end;
                end;
                lastMol := rMax[k+256];
                total := 0;
            end;
            i := i + 1;
            rMin[i] := k;
            rMax[i] := rMax[k+256];
            rArea[i] := area;
            rUser1[i] := sum;
            rAngle[i] := sum - area*AvgBkg;
            total := total + rAngle[i];
            rX[i] := rX[k+256];
            rY[i] := rY[k+256];
        end;
    end;
    {SetCounter(i);} {zeros measurements above i}
    disposePic(sumPid);
    SetUser1Label('sum');
    SetUser2Label('s-a*b/t');
    ShowResults;
    RestoreState;
    SelectPic(Seg8bPid);
    ShowResults;
    UpdateResults;
    Copy; {manually delete the junk at end}
end;

macro 'make stack';
begin
    setNewSize(512,100);
    makeNewStack('stack');
    stackPid := pidNumber;
    setMemo('stackPid',stackPid);
end;

macro 'front image is stack';
begin
    stackPid := pidNumber;
    setMemo('stackPid',stackPid);
end;

macro '[F15]copy molecule to stack';
var
    left, top, width, height: integer;
    tempPid, srPid: integer;
    pwidth, pheight: integer;
    swidth, sheight: integer;
    s: string;
begin
    SetCounter(513);
    SetOptions('X-Y Center,angle');
    Measure;
    if rAngle[rCount] > 90 then rAngle[rCount] := rAngle[rCount] - 180;
    InsetRoi(-10);
    GetRoi(left, top, width, height);
    if width < height then begin
        top := top - height/4;
        height := height + height/2;
    end else begin
        left := left - width / 4;
        width := width + width/2;
    end;
    MakeRoi(left, top,width, height);
    InsetRoi(-10);
    Duplicate('temp');
    tempPid := pidNumber;
    if width < height then begin
        rotateLeft(true);
        rAngle[rCount] := rAngle[rCount] + 90;
        srPid := pidNumber;
        if pidExists(tempPid) then begin
            choosePic(tempPid);
            dispose;
        end;
        tempPid := srPid;
        choosePic(tempPid);
    end;
    SetScaling('Bilinear,New Window');
    {SetScaling('New Window');}
    ScaleAndRotate(1,1,rAngle[rCount]);
    srPid := pidNumber;
    ChoosePic(tempPid);
    Dispose;
    SelectPic(srPid);
    GetPicSize(pWIdth, pHeight);
    ChoosePic(stackPid);
    GetPicSize(sWidth, sHeight);
    SelectPic(srPid);
    SelectAll;
    Copy;
    Dispose;
    SelectPic(stackPid);
    SelectSlice(nSlices);
    AddSlice;
    Paste;
    ChoosePic(raw16Pid);
    s := GetPicName;
    ChoosePic(stackPid);
    MoveTo(10,10);
    WriteLn(s);
    WriteLn(top);
    WriteLn(left);
    SelectPic(Proc8Pid);
    SetDensitySlice(255,255);
end;